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1.  Abstract 


Objective: 

The  objective  of  this  proposal  is  to  develop  physically  based  models  to  predict  the  penetration 
depth  of  common  military  munitions  in  various  soil  conditions.  Ultimately,  the  models  will  be 
used  to  determine  probable  depths  of  munitions  in  the  soil  of  formerly  used  defense  sites  in  support 
of  planning  for  remediation.  The  simulation  results  can  be  used  to  aid  sensor  detection  and  removal 
of  these  munitions. 

Technical  Approach: 

To  model  munitions  penetration,  a  meshfree  framework  based  on  the  Reproducing  Kernel  Particle 
Method  (RKPM)  is  developed  for  handling  extremely  large  deformation.  A  particle-to-particle 
based  contact  algorithm  is  introduced  to  efficiently  capture  energy  and  momentum  exchanges 
between  munitions  and  soils  and  between  soil  fragments.  A  two-field  (displacement  and  pressure) 
semi-Lagrangian  formulation  is  developed  and  implemented  in  consideration  of  porous  nature  of 
soils.  Stabilized  nodal  domain  integration  schemes  that  ensure  the  accuracy  and  stability  of  the 
numerical  solution  is  developed  for  the  two-field  Galerkin  formulation. 

To  accurately  represent  the  behavior  of  the  soil,  a  viscoplasticity  model  is  developed  with 
regularized  softening  to  account  for  large  deformation  of  the  soils.  The  model  accounts  for  the 
behaviors  seen  in  penetrations  problems,  including  nonlinear  pressure  sensitivity  of  shear  strength, 
rate  dependence,  shear-enhanced  dilation,  and  a  compaction  hardening  due  to  pore  collapse  and 
grain  crushing.  The  model  is  updated  with  a  regularized  softening  with  increasing  porosity  that 
can  naturally  transition  to  more  fluid-like  flow,  including  liquefaction  that  is  sometimes  observed 
during  penetration.  The  viscoplasticity  model  is  embedded  in  a  saturated  two-field  meshfree  code, 
with  partially  saturated  framework  to  be  completed  in  subsequent  research. 

Results: 

The  two-field  (displacement-pressure)  formulation  based  on  the  Biot  theory  has  been  developed  and 
implemented  under  the  semi-Lagrangian  RK  framework,  where  displacement  and  pressure  field  are 
independently  approximated  by  the  semi-Lagrangian  Reproducing  Kernel  (RK)  shape  functions. 
Some  numerical  schemes  originally  designed  for  the  single-field  formulation  have  been  modified  and 
implemented  for  the  two-field  formulation,  including  the  modified  stabilized  non-conforming  nodal 
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integration  for  the  domain  integration,  stress  update,  and  kernel  contact  algorithms.  The  central 
difference  and  forward  Euler  temporal  integration  schemes  have  been  applied  to  the  displacement 
and  pressure  fields,  respectively,  in  the  two-field  formulation,  leading  to  an  explicit  time  marching 
scheme. 

The  three-invariant  viscoplasticity  model  has  been  developed  to  effectively  integrate  tensile, 
shear,  and  compressive  behavior.  The  evolution  of  volumetric  plastic  strain  has  been  explicitly 
connected  to  the  void  ratio,  allowing  the  model  to  be  integrated  with  a  poromechanical 
framework.  Novel  hardening/softening  laws  have  been  added  to  characterize  strengthening  and 
weakening  in  different  loading  regimes,  and  regularized  the  softening  using  viscoplasticity.  The 
model  also  accounts  for  rate  effects,  differences  in  strength  in  triaxial  extension  and 
compression,  compression  hardening,  and  other  effects.  An  efficient  implementation  using  the 
spectral  decomposition  has  been  employed  to  reduce  the  cost  of  the  complicated  model. 

The  model  has  been  verified  against  a  Drucker-Prager  model,  and  the  meshfree  and  finite 
element  implementations  have  been  verified  against  each  other  to  ensure  proper  implementation. 
The  behaviors  of  the  model  have  been  demonstrated  in  a  numerical  framework  using  reasonably 
simple  example  problems.  The  developed  two-field  meshfree  code  have  been  employed  to  simulate 
penetration  process  into  soil  and  predict  the  final  penetration  depth  under  different  penetration 
angles.  In  the  penetration  simulations,  the  numerical  results  show  that  the  maximum  penetration 
depth  varies  from  3m  to  lm  with  penetration  angles  ranging  from  90°  to  45°.  The  deformation  of  soil, 
e.g.  soil  splashing  on  the  free  surface,  reflects  the  experiment  observations  after  impact. 

Benefits: 

The  project  will  benefit  Department  of  Defense  by  providing  a  robust  numerical  framework  for 
modeling  penetration  into  soils.  Eventually,  the  results  of  the  simulations  will  translate  into  a  set 
of  tables  for  probable  depths  of  munitions  based  on  soil  conditions,  projectile  type,  and  firing 
conditions. 

The  modeling  will  also  have  broader  impacts  to  the  engineering  and  scientific  communities. 
The  framework  will  be  able  to  model  other  penetration  scenarios  for  soil,  rock,  and  concrete,  for 
applications  as  diverse  as  deep  penetrators  designed  to  target  underground  bunkers  to  meteor 
impacts  on  extraterrestrial  bodies.  The  constitutive  models  will  further  our  fundamental 
understanding  of  soil  behavior  as  we  move  toward  physically  based  models  for  capturing  observed 
soil  responses.  In  addition,  the  numerical  algorithms  will  enhance  the  set  of  tools  available  to  solve 
many  physical  problems,  especially  those  involved  large  deformation,  material  separation,  and 
coupled  physics  problems. 
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2.  Objective 

The  objective  of  this  project  is  aimed  to  develop  a  meshfree  method  for  soil-munitions 
interaction  in  order  to  predict  penetration  depth  for  a  range  of  characteristics.  The  numerical 
framework  for  modeling  the  impact/penetration  is  based  on  the  RKPM.  This  method  has  the 
advantages  that  it  can  model  extreme  distortion  and  material  separation  with  relative  ease  and 
that  it  can  directly  incorporate  the  physically  based  soil  constitutive  models.  Moreover,  as  a 
continuum  approach,  the  computation  cost  in  the  current  framework  is  more  manageable  in 
comparison  to  discrete  element  methods.  These  features  are  essential  for  creating  an  accurate, 
tractable,  and  physically  based  model  of  projectile  penetration. 

Equally  important  to  developing  a  broadly  applicable  computing  framework  is  modeling  the 
material  behavior.  A  constitutive  model  for  the  soil,  a  three-invariant  viscoplasticity  model,  is 
adapted  to  reproduce  the  physical  behaviors  of  soils,  including  hardening,  softening,  and  rate 
dependency.  This  solid  soil  model  is  currently  implemented  in  the  displacement  based 
formulation,  and  will  be  embedded  in  a  poromechanical  framework  to  reproduce  the  effects  of 
water  in  the  soil  in  the  future. 

Ultimately,  the  developed  meshfree  framework  with  the  soil  model  is  used  to  simulate 
penetration  of  munitions  into  soils  for  a  range  of  soil  and  munitions  parameters,  velocities  and 
impact  angles.  In  the  future,  the  outputs  of  these  simulations  will  be  used  create  tables  or 
regression  curves  of  probable  depths  based  on  soil  parameters  that  can  be  measured  by  standard 
laboratory  tests,  known  munitions  characteristics,  and  the  scenarios  under  which  those  projectiles 
were  launched. 

The  soil  model  is  developed  from  the  Sandia  Geomodel  (Fossum  and  Brannon  2004a  and 
2004b,  Foster  et  al.,  2005,  Motamedi  and  Foster,  2017)  This  model  that  accounts  for  many 
aspects  of  soil  mechanical  behavior,  including  rate  dependence,  difference  in  triaxial  extension 
and  compression  strength,  kinematic  hardening  and  nonassociative  dilation  at  low  mean  stresses, 
compaction  hardening  at  high  mean  stress,  and  other  features.  The  model  will  need  to  be 
extended  to  include  poromechanical  effects  under  saturated  conditions  (unsaturated  behavior  will 
be  left  for  future  work)  directly  related  to  material  porosity.  The  model  will  also  be  extended  to 
the  finite  deformation  regime. 

The  completed  model  will  be  validated  against  both  laboratory  experiments  (Seguin,  et  al. 
2008)  and  field  data  (Christiensen,  1989)  both  for  dry  and  saturated  soils.  The  models  will  be 
adjusted  until  reasonable  predictive  agreement  about  penetration  depths,  within  20%  relative 
error,  can  be  obtained  over  a  variety  of  conditions.  These  models  will  serve  as  a  baseline  for  a 
future  extension  to  partial  saturation,  and  eventually  the  production  of  tables  for  determination  of 
depth  for  a  variety  of  conditions. 
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3.  Background 

The  penetration  of  projectiles  into  the  earth  is  a  widely  studied  but  complex  problem.  The  impact 
of  bodies  into  soil  and  rock  has  a  wide  variety  of  applications,  from  the  study  of  meteor  impacts 
to  military  applications  to  civilian  engineering  problems  such  as  deep  compaction.  In  military 
applications,  there  are  a  wide  variety  of  scenarios,  from  deep  penetrators  designed  for 
underground  bunkers  to  impacts  of  accidental  crashes. 

In  this  objective  of  the  SEED  project  is  to  examine  the  penetration  of  munitions  at  formerly 
used  military  test  sites.  The  purpose  of  this  examination  is  to  determine  the  probable  depth  of 
buried  projectiles  so  that  they  may  be  easily  detected  and  removed  as  the  site  is  remediated.  The 
determination  of  the  final  depth  is  a  complex  task.  It  depends  soil  constituents:  grains  of  different 
sizes,  shapes,  chemical  cohesion,  and  mineral  content,  as  well  as  the  soil  state,  including  porosity 
and  degree  of  saturation.  Some  of  these  properties  change  locally  during  an  impact  event  through 
compaction,  shearing,  and  grain  crushing.  The  depth  also  depends  on  the  shape,  mass,  and 
mechanical  properties  of  the  projectiles.  Finally,  the  impact  angle,  velocity  and  angular  velocity 
greatly  influence  the  ultimate  depth.  (Omidvar  et  al.,  2014) 

Even  well  characterized,  the  ultimate  outcome  of  a  penetration  event  can  be  difficult  to 
predict.  The  physics  are  complicated,  and  no  analytical  method  could  hope  to  accurately  predict 
the  outcome.  Numerical  techniques  may  be  employed,  but  extreme  deformation  and  complex 
physics  makes  the  problem  extremely  challenging. 

Experimental  studies  of  penetration  can  be  divided  two  categories:  laboratory  and  field 
experiments.  While  laboratory  experiments  are  more  controlled  and  more  data  can  be  readily 
measured,  they  are  generally  limited  in  scale,  and  can  be  difficult  to  extrapolate  to  field-scale 
simulations.  Despite  greater  uncertainty  in  material  state,  velocity,  and  other  parameters  in  field 
experiments  (which  have  similar  uncertainty  to  sites  to  be  remediated),  they  remain  vital  to 
validation.  Laboratory  experiments  of  penetration  are  subject  to  confining  boundary  conditions 
that  can  greatly  affect  penetration  depth  and  the  type  of  deformation  of  the  soil  particles  (Seguin, 
et  al  2008,  Omidvar,  et  al  2014).  The  interactions  are  complex;  however,  and  sometimes  give 
difficult  to  interpret  results  (Borg,  et  al  2013).  Ideally,  fully  transmitting  boundary  conditions 
would  allow  stress  waves  to  propagate  away  from  the  projectile  without  reflection.  Such 
boundary  conditions  are  difficult  to  set  up  experimentally,  however.  Fortunately,  numerical 
models  can  account  for  varying  boundary  conditions,  and  a  physically  based  model  can  be 
validated  in  a  controlled  setting  before  being  applied  to  more  open  field  experiments.  An 
excellent  compilation  of  field  data  is  the  SAND  Report,  Twenty  five  years  of  penetration  records 
at  Sandia  National  Labor atories-PENTDB:  a  relational  database,  (Christiansen,  1989)  though 
this  document  is  limited  release. 

While,  experiments  are  vital  to  understanding  the  behavior  of  penetration  in  soils,  it  is 
impossible  to  run  experiments  for  every  scenario  at  a  test.  Modeling  is  a  cost-effective  solution 
to  capturing  the  response  of  projectile  impact  in  soils.  With  regards  to  penetration  depth,  a 
number  of  phenomenological  and  semi-empirical  methods  have  been  developed,  e.g.  (Euler, 
1914;  Carter  et  al.,  1986;  Forrestal  and  Luk,  1992,  Durbin  and  Masri,  2004;  Salgado,  et  al.,  1997; 
Yu  and  Houlsby,  1991).  Many  of  these  models  are  reviewed  in  Omidvar,  et  al.,  2014.  These 
models  vary  in  the  extent  of  the  physics  involved  as  well  as  the  extent  of  validation,  but  it  is 
impossible  to  solve  all  of  the  equations  physically  in  a  very  general  way  that  can  be  easily 
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extended  to  a  variety  of  situations.  Still,  the  models  can  provide  good  initial  verification  to 
numerical  codes. 

The  equations  for  physics-based  modeling  of  penetration  are  far  too  complex  to  be  solved  by 
analytical  methods.  Some  simulations  have  been  performed  with  finite  element  analysis.  While 
the  finite  element  method  (FEM)  has  been  used  for  a  wide  variety  of  mechanical  and  structural 
problems  with  great  success,  the  penetration  problem  is  not  as  easily  solved.  Due  to  challenges 
with  material  separation  a  narrow  initial  cavity  is  often  assumed.  Other  challenges  include  very 
large  distortion  of  the  elements,  which  can  lead  to  inaccurate  results,  and  the  evolving  contact 
problem.  Though  a  number  of  algorithms  have  been  developed  for  contact,  practical  impact 
problems  remain  a  challenge. 

Recently,  particle-based  models,  notably  discrete  or  distinct  element  methods  (DEM),  have 
been  used  to  examine  soil  deformation  in  the  case  of  projectiles.  Here  the  particles  are  modeled 
individually,  usually  with  simple  constitutive  laws  governing  the  stiffness  between  the  individual 
particles.  While  computationally  expensive,  they  have  been  used  to  model  penetration  in  sands 
and  gravels,  e.g.  (Tiwari  et  al.,  2015;  Borvik  et  al.,  2015).  Some  approximate  scaling  laws  are 
required  for  fine-grained  soils,  though  it  remains  unclear  how  generally  these  apply  to  a  range  of 
loading  scenarios.  While  the  approach  can  capture  some  mesoscale  phenomena  such  as  the 
formation  and  buckling  of  force  chains,  challenges  such  as  physically  based  grain  fracture  and 
sophisticated  frictional  interaction  and  wear  of  real  grains  remain.  DEM  have  often  been  coupled 
with  rigid  or  finite  element  models  of  the  penetrators,  e.g.  Onate  and  Rojek,  2004. 

As  opposed  to  DEM,  meshfree  methods  based  on  the  continuum  theory,  such  as  the  Element- 
Free  Galerkin  method  (Belytschko  et  al.,  1994)  and  the  RKPM  (Liu  et  al.,  1995)  have  gained 
popularity  in  modeling  large  deformation  problems.  Chen  et  al.  (1996)  introduced  a  consistent 
discrete  formulation  for  RK  approximation  and  extended  the  RK  formulation  to  finite 
deformation  problems.  Guan  et  al.  (2014)  introduced  the  semi-Lagrangian  RKPM  formulation 
with  the  kernel  contact  algorithm  under  the  meshfree  framework.  Bessa  et  al.  (2014)  established 
a  direct  link  between  the  RKPM  and  the  state-based  peridynamics  methods  (Silling,  2000;  Silling 
&  Lehoucq,  2008)  and  concluded  that  the  state-based  peridynamics  leads  to  an  approximation  of 
the  derivatives  that  can  be  obtained  from  the  RKPM.  In  the  Galerkin  numerical  procedure, 
meshfree  methods  require  high-order  numerical  integration  on  a  background  grid,  resulting  in 
considerable  computation  efforts  (Dolbow  and  Belytschko,  1999).  On  the  other  hand,  low-order 
numerical  integration,  such  as  nodal  integration,  leads  to  instability.  In  order  to  circumvent  the 
numerical  instability,  Bonet  et  al.  (1999)  introduced  a  correction  term  into  the  derivative  of  shape 
function  at  nodal  point  by  enforcing  the  linear  exactness  condition.  Chen  et  al.  (2001)  proposed 
the  Stabilized  Conforming  Nodal  Integration  (SCNI)  in  which  the  linear  exactness  is  achieved  in 
the  Galerkin  approximation  of  second  order  PDEs,  and  later  applied  SCNI  to  large  deformation 
problems  in  Chen  et  al.  (2002).  Recently,  Chen  et  al.  (2013)  generalized  the  SCNI  and  the  idea 
of  integration  constraint  and  developed  a  variationally  consistent  integration  scheme  to  achieve 
arbitrary  order  of  exactness  in  the  Galerkin  formulation.  A  non-confirming  counterpart  of  the 
SCNI,  where  the  smoothed  strain  is  constructed  over  non-confirming  subdomains,  has  been 
introduced  in  concrete  penetration  modeling  (Guan  et  al.,  2011;  Chi  et  al.,  2015).  Hillman  et  al. 
(2014)  introduced  stabilized  integration  schemes  under  the  variationally  consistent  integration 
(VCI)  framework  for  modeling  impact  and  penetration  problems.  An  implicit  gradient 
stabilization  was  later  introduced  for  enhanced  stability  in  the  Galerkin  meshfree  method  when 


5 


the  nodal  integration  is  used  (Hillman  and  Chen,  2016).  Yreux  and  Chen  (2017)  introduced  a 
new  formulation  to  construct  the  Moving  Least  Squares\RK  approximation  even  when  the 
neighboring  nodes  are  insufficient  due  to  material  separation. 

The  history  of  modern  soil  mechanics  is  relatively  young,  reaching  back  only  a  hundred 
years  or  so  to  Terzhaghi’s  effective  stress  principle  and  one-dimensional  consolidation  theory.  In 
the  latter  half  of  the  twentieth  century  and  into  the  twenty-first,  however,  and  profusion  of 
models  have  been  introduced  to  capture  various  behaviors  of  soils. 

Among  continuum  models,  plasticity  models  and  their  extensions  have  been  applied  to  a 
variety  of  applications  with  success.  One  of  the  first  models  was  the  Mohr-Coulomb  model, 
developed  from  a  failure  condition  for  soils,  and  captures  the  mean  stress-dependent  shear 
strength  of  soil.  Others,  such  as  the  Drucker-Prager  and  Lade-Duncan  models,  built  on  these 
characteristics  while  improving  computational  efficiency.  The  Drucker-Prager  criterion  (Drucker 
and  Prager  1952),  among  others,  is  a  suitable  candidate  as  a  preliminary  constitutive  model  for 
many  engineering  applications  due  to  its  simplicity,  smooth  yield  surface,  relatively  easy 
implementation  and  high  efficiency  within  numerical  algorithms.  However,  it  has  a  few 
limitations,  including  overestimating  strength  in  triaxial  extension  (Alejano,  2012)  and  being 
unable  to  model  soil  compaction  at  high  mean  stresses. 

The  original  and  modified  Cam-Clay  models  (Roscoe  and  Burland,  1968),  do  account  for 
both  dilation  and  compaction  and  different  mean  stresses,  but  have  limitations  in  the  shape  of  the 
yield  surface  and  are  designed  for  slow  consolidation  behavior.  A  family  of  cap  models  has  been 
developed  for  more  general  response  of  both  soils  and  rock,  some  of  which  are  discussed  in 
(Fossum  and  Brannon  2004a).  An  enormous  number  of  models  have  been  developed  to  simulate 
various  soil  behaviors.  An  early  review  of  soil  plasticity  models  may  be  found  in  (Scott  1985), 
among  other  places.  Important  characteristics  of  accurate  models  include  pressure  dependence, 
which  is  nonlinear  over  large  ranges,  differences  in  triaxial  extension  and  compression  strength, 
and,  at  large  pressures,  the  ability  to  capture  inelastic  compaction  due  to  pore  collapse  and  grain 
fracture.  Rate  dependent  effects  become  important  for  applications  involving  high  strain  rates, 
including  penetration. 
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4.  Materials  and  Methods 


4.1.  Semi-Lagrangian  RK  framework  for  penetration  modeling 

4.1.1.  Semi-Lagrangian  RK  Approximation 

Penetration  simulations  involve  extremely  large  deformation  and  highly  fragmented 
configurations  that  are  very  difficult  to  capture  with  conventional  finite  element  methods. 
Therefore,  a  meshfree  framework  based  on  the  semi-Lagrangian  RKPM  (Guan  et  al.,  201 1;  Chi  et 
al.,  2015,  Chen  et  al.,  2017)  is  introduced  to  effectively  model  excessive  soil  deformation. 

Consider  a  set  of  NP  scattered  points  in  a  closed  domain,  Q.x  =  Cix  u  dQx ,  with  spatial  positions 

in  the  initial  configuration  where  is  the  boundary  of  the  open  domain  in  the 

initial  configuration.  In  this  report,  X  and  X  are  designated  for  coordinates  in  the  initial 
configuration  whereas  x  and  x  for  coordinates  in  the  current/deformed  configuration.  The  two 
coordinates  are  related  through  a  mapping  x  =  <p(X,t),  and  the  material  point  in  the  current 

configuration  is  defined  as  x7  =  <p(Xnt  ) .  The  semi-Lagrangian  RK  shape  function  is  constructed 
following  a  material  point  xi  while  its  domain  of  influence  has  fixed  size  and  shape  (Figure  1).  The 
semi-Lagrangian  RK  approximation  of  a  function  u(x)  is  expressed  as: 

NP  NP 

=Zc(x;x-x/K(x-x/K  •  (!) 

7=1  7=1 

where  VF/  is  the  semi-Lagrangian  shape  function  associated  with  node  I  and  di  is  the 
corresponding  coefficient.  Above  the  kernel  function  (j>a  (x-x,)  controls  the  smoothness  and  the 

domain  of  influence  of  the  shape  function  with  a  support  size  a.  The  center  of  the  kernel  function 
follows  the  material  point  while  its  domain  of  influence  is  fixed  in  space.  The  correction  function 
C(x;x-x7)  is  introduced  to  ensure  solution  accuracy  up  to  the  desired  order,  and  it  is  usually 
expressed  by  the  linear  combination  of  a  set  of  m-th  order  complete  monomials,  that  is 

m 

C(x;x-x/)=  ^(x-x/)/?^(x)  =  H(x-x/)b(x)  (2) 

1/71=0 

where  /?  is  the  multi-dimensional  index,  considering  up  to  three  dimensions:  Jd  =  (j3vj32,j3i) 
with  its  length  defined  as  \p\  =  ,  bfj  =  is  the  corresponding  coefficient 

(x-x^  =  (Xj  -xuYl{x2  -x2/)^2(x3  -x3IY}  ,  and  b(x)  and  H(x-x/)  are  the  vector  forms  of 
bp  and  (x-x7)^. 

The  coefficient  b(x)  can  be  determined  by  the  reproducing  condition: 

NP 

^T//(x)xf  =x/s  ;  \p\<m,  (3) 

/= i 
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and  subsequently  the  RK  shape  function  is  given  as  follows. 

(x)  =  Hr  (O)lVr 1  (x)H  (x  -  x, )  (f)a  (x  -  x, )  (4) 

where 

NP 

M(x)  =  Yj  h (x - x,  )HT  (x - x,  )</>a  (x - x, )  (5) 

/=i 

It  is  noted  that  the  kernel  function  abovementioned  follows  the  material  points  while  its  domain 
of  influence  is  fixed,  which  has  been  called  the  Eulerian  kernel  in  literature.  The  influence  of  the 
Eulerian  kernel  can  serve  for  contact  detection  in  the  points-based  simulations  (Li  et  al.,  2001; 
Guan  et  al.,  2011,  Chi  et  al.,  2015),  and  the  semi-Lagrangian  RK  shape  function  inherits  such 
property.  Thus,  several  contacts  algorithms  have  been  introduced  to  impose  contact  conditions 
without  a  pre-defined  contact  surfaces  under  the  semi-Lagrangian  RK  framework  (Chi  et  al, 

2015;  Guan  et  al.,  2011).  Specially,  a  levelset-enhanced  frictional  kernel  contact  algorithm  is 
adopted  in  this  project.  The  details  are  referred  to  (Chi  et  al.,  2015). 
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Figure  1.  Semi-Lagrangian  Reproducing  Kernel  (RK)  and  meshfree  discretization 

4.1.2.  Single-field  Semi-Lagrangian  RK  Formulation 

The  equation  of  motion  of  a  deformable  body  in  the  current/deformed  configuration  with  traction 
and  displacement  boundary  conditions  is  given: 


pix  =  V  o  +  b 

in  Qx 

no  =  h 

on  dQhx 

(6) 

u  =  g 

on  dQ.gx 

where  p  is  the  density;  b  is  the  body  force;  u  is  the  displacement  vector;  o  is  the  Cauchy  stress; 

is  the  open  domain  in  the  current  configuration  with  its  enclosed  boundary  dClx ;  dQx  and 
dQx  are  the  boundaries  with  the  prescribed  traction  h  and  the  prescribed  displacement  g, 
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respectively;  and  dQx  u  dQx  =  dQx ,  dQx  n  dQx  =  0 .  The  variational  equation  of  equation  (6) 
is: 


f  pdu-iidQ  +  f  SVu  ■  add  =  f  <5u  b<iQ+f  ,  <5u-h<TT.  (7) 

Jn*  Jqx  Jqx  Jdahx 

The  displacement  u  and  its  variation  £>tl  are  approximated  using  the  semi-Lagrangian  RK  as  in 
equation  (1),  and  the  associated  velocity  and  acceleration  are  approximated  by  the  material  time 
derivatives  of  equation  (1),  leading  to: 


u"  (X,  f )  =  V*  (X,  f)  =  £  (^  (x)d(O  +  V,  (x)d(O) , 

7=1 

ii'  (x,  t )  =  ah (x,  t )  =  £('P/  (x)d(t)  +  2V,  (x)d(f)  +  f  7  (x)d(f))  • 


(8) 


(9) 


Above  the  superscript  h  denotes  the  finite  dimensional  approximation;  d(r) ,  d(r) ,  and  d(t)  are 

the  nodal  coefficients  of  displacement,  velocity,  and  acceleration,  respectively;  'F  and  VP  are 
the  time  rate  changes  due  to  the  Eulerian  kernel  function,  i.e., 


<4(x-x/)=^ 


x-x. 


((X_X/).(V_V/) 


a  x-x. 


(10) 


Substituting  equations  (8)  and  (9)  into  the  Galerkin  form  of  equation  (7),  it  leads  to  the  semi¬ 
discrete  equation  in  the  matrix  form  as  follows. 


where 


Md  (t )  +  Gd  (t)  +  Gd  (t)  =  Fext  -  Fint , 

(ii) 

M„=Jn  /W.Vjidn, 

(12) 

G„=Jn  2{H’pJldn, 

(13) 

(14) 

F;-v'  =  [  T.bJQ+f  V,hdr, 

1  Jci;  1  Jen1’  1 

(15) 

f;iu  =  b'  Ec/q  . 

(16) 

Above  I  denotes  the  identity  matrix,  B,  is  the  gradient  matrix  representing  the  symmetric  part  of 

Vu  associated  with  node  /,  L  is  the  stress  vector  associated  with  o .  The  matrices  associated 
with  the  time  rate  change  of  the  semi-Lagrangian  kernels,  G ,  G ,  are  usually  omitted  for 
computational  efficiency  since  they  have  very  subtle  effects  on  the  solutions  when  the  nodal 
integration  is  used  (Chi  et  al.,  2015).  Numerical  time  integration  schemes,  such  as  the  Newmark 
methods,  can  then  be  introduced  in  equation  (11)  (see  Siriaksorn  et  al.,  2017  for  details).  In  this 
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project  the  central  difference  is  used  for  temporal  integration  for  the  displacement  unless 
otherwise  noted. 


4.1.3.  Two-field  (u-p)  Semi-Lagrangian  RK  Formulation 


To  account  for  the  coupling  effect  between  solid  and  fluid,  a  two-field  meshfree  formulation 
based  on  the  Biot  mixture  theory  (Biot,  1956)  is  developed  and  implemented.  At  the  current 
stage,  only  fully  saturated  soil  material  is  considered.  The  governing  equations  are  the  balance  of 
momentum  and  continuity  equation  for  the  fluid  phase,  respectively: 

/?ii=Vo+b  , 

in  Qx 

(17) 

„  .  P  „ 

aV -uh - i-V  -q 

M 

=  0 ,  in  Qx 

(18) 

with  boundary  conditions: 

no  =  h 

on  oOx 

u  =  g 

-n  q  =  v. 

on  dQf 
on  dQx 

(19) 

P  =  P 

on  sq; 

where  p  is  the  saturated  density  of  the  porous  medium,  a  is  the  Biot  coefficient,  P  is  the  pore 

pressure,  M  is  the  Biot  compressibility  modulus,  and  q  is  the  superficial  velocity,  vs  is  the 
prescribed  fluid  inflow  on  dQx ,  P  is  the  prescribed  pore  fluid  pressure  on  oQ' ,  and 


dQx  u  dQx  =  dQx ,  dQ.sx  n  BQ'x  =  0  .  The  total  stress  can  be  decomposed  into  partial  stresses 
from  solid  and  fluid,  i.e.,  o  =  (1  -  X)os  +  XcrFI ,  where  as  is  the  solid  stress,  crF  is  the  fluid 
stress,  X  is  the  porosity,  defined  as  the  ratio  of  the  void  volume  to  the  total  volume  V.  Introducing 
the  Biot  coefficient  a,  the  total  stress  can  be  alternatively  represented  as 


o  =  o  -  aPl ,  (20) 

where  o  is  the  effective  stress  of  the  solid  phase.  The  soil  constitutive  law  described  in  the 
following  section  will  be  introduced  for  determining  the  effective  stress  o  .  The  Biot  coefficient 
a  and  Biot  compressibility  modulus  M  can  be  defined  respectively  as 


<2  =  1- 


K 


1  _  a-X  X 
M  ~  Ks  K/ 


(21) 

(22) 


where  K  is  the  bulk  modulus  of  the  porous  medium;  Ks  is  the  bulk  modulus  of  the  solid  grains; 
Kf  is  the  bulk  modulus  of  the  fluid.  For  isotropic  porous  media,  the  relationship  between  the  flux 
or  superficial  velocity  q  and  the  pore  pressure  P  can  be  described  by  using  Darcy’s  law  as 
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q=-— (vp-pfg), 

jUx  ’ 


(23) 


where  k  is  the  intrinsic  permeability,  p  is  the  fluid  dynamic  viscosity,  Pf  is  the  fluid  density,  and 
g  is  the  gravity. 

The  variational  form,  and  subsequently  the  Galerkin  formulation,  of  (17)  and  (18)  are  given 
as  follows. 


f  SVuh:ahdCl-{  SV  uhaPhdQ.  +  f  5nh  ■  pnh  dQ.=  \  ,8nh  hdT +[  5uh -bdCL,  (24) 

JC2X  JQX  JQX  JdQx  Jax 


f  8PhdV  uhdn-\  8Ph  —dQ.+  f  —SVPh -VPh  dfi. 

J nx  Jnx  m  Jnx  ^ 

=  f  SP\dT+\  —pf8VPh  gdCl 

Jan*  5  Jnx 

(25) 

The  central  difference  and  the  forward  Euler  temporal  integration  schemes  are  employed  for  the 
displacement  field  and  the  pressure  field,  respectively,  as  follows. 

d;+1  =d';+Atv;+0.5At2a;, 

(26) 

v”+1  =  v"+1  +  0.5Am"+1 , 

(27) 

Pi  ~  Pi  ^tpi  » 

(28) 

where  At  is  the  time  step,  the  superscript  n  denotes  the  time  step  count,  d,  v,  a,  and  p  are  the  nodal 
coefficients  for  displacement,  velocity,  acceleration,  and  fluid  pressure,  respectively.  The 
predicted  velocity  v"+1  is  defined  as 

v?+1  =v;+0.5A /a” 

(29) 

After  introducing  the  semi-Lagrangian  RK  shape  function  (equation  (1))  for  the  approximation 
functions  of  displacements  uh  and  fluid  pressure  Ph  and  their  time  derivatives  using  equations 
(8)  and  (9),  the  equations  (24)  and  (25)  are  recast  into  the  matrix  forms: 

Ma"+1  +  Gv',+1  +  Gd',+1  =  Fext  -  Fint , 

(30) 

Sp"+1  +  Qp"+1  +  Qd"+1  =  Fext  -  Fint , 

(31) 

where  M,  G,  G ,  and  Fext  are  the  same  as  in  equations  (12)-(15),  and 

(32) 

(33) 

Q„=ln  aV'VV.dn, 

(34) 
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F;nt  =  J  VT,  •  Gn+1dQ.  -  J  ccV'VjP^dn  -  J  T,brfQ ,  (35) 

f;xi  =  f  T>^r,  (36) 

J  0QX  1  s 

F/nt  =  f  V  •  u"+1  dn  +  J  (k  /  ju) VV,  ■  VPn+1dn  -\(k/p)  pf' VV ,  ■  gdQ. .  (37) 

n 

Applying  the  central  difference  and  forward  Euler  schemes  in  equations  (26)-(29)  leads  to  the 
full  discrete  equations: 

(M  +  0.5Af(C  +  G))a"+1  =  Fext  -  Fint  -  (C  +  G)vn+1  -  Gd"+1 ,  (38) 

Sp"+1  =  Fext  -  Fint  -  Qd"+1  -  Qp"+1 .  (39) 


Above,  a  mass  proportional  damping  matrix  C  can  be  introduced  if  desired.  Similar  to  Section 
4.2,  the  matrices  associated  with  the  time  rate  change  of  the  semi-Lagrangian  kernels,  G ,  G , 

G ,  and  G ,  can  be  omitted  for  computational  efficiency.  Subsequently,  equations  (38)  and  (39) 
degenerate  to 


(M  +  0.5AtC)a"+1  =  Fext  -Fint  -Cvn+1 , 

(40) 

Sp”+1  =  Fext  -  Fint . 

(41) 

The  lumped  matrix  scheme  by  the  row  sum  is  employed  for  M  and  S  to  acquire  diagonal 
matrices.  This  process  leads  to  an  explicit  scheme,  the  information  of  at  the  (n+1  )-th  time  step  is 
fully  independent  on  the  information  from  the  n-th  time  step,  and  therefore  it  is  efficient  for 
highly  dynamic  applications.  However,  the  scheme  is  conditionally  stable.  The  stable  time  step 
size  can  be  estimated  using  the  von  Neumann  method  and  the  details  are  referred  to  (Chi  and 
Siriaksom,  2017). 


4. 1 .4.  Domain  integration  in  the  Galerkin  formulation 

The  domain  integration  in  the  meshfree  Galerkin  formulation  (equations  (7)  and  (24)-(25)) 
requires  special  attention  due  to  the  nature  of  the  RK  approximation  and  the  applications  involving 
extreme  deformation  and  discontinuities  (Chen  et  al.  2013;  Hillman  et  al.,  2014;  Chen  et  al.  2017). 
In  this  project,  the  Modified  Stabilized  Nonconforming  Nodal  Integration  (MSNNI)  is 
implemented  (Chen  et  al.,  2006).  The  direct  gradients  involved  in  equations  (12)-(16)  and  (32)- 
(37)  are  replaced  by  an  assumed  gradient  operator  to  construct  smoothed  derivatives  of  the  shape 
function  over  a  prismatic  nodal  representative  domain  (Figure  2),  that  is 

V'P,(xI)  =  -f|r  (x)n(x)rfr  (42) 

L 

where  VL  and  T7  are  the  volume  and  boundary,  respectively,  of  the  nodal  representative  domain 
of  node  L  and  n  is  the  unit  outward  normal  vector  of  T7 .  Therefore,  the  smoothed  strain  field 
and  pressure  gradient  field  can  be  constructed  as: 
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(43) 


NP  Nsl  _  _  t  —  — 

wZZ(B/(x^-B/(xp)  E(By(xJ-B/xPK  (48) 

L= 1  £=1 

NP  NsL  _  _  i  _  _ 

-(B!(xJ-B'(x{))v{  (49) 

L=1  i'=l  H 

are  added  in  M  and  S,  respectively,  in  equations  (40)  and  (41).  Above  w  and  wp  are  stabilization 
parameters  ranging  between  0  and  1,  E  is  the  elastic  part  of  the  material  tangent  tensor,  Vc  is  the 

nodal  volume  associated  with  subdomain  j  (see  Figure  2),  V)  =  VL ,  and  NsL  is  the  number 

of  subdomains  in  the  nodal  representative  domain  of  node  L.  The  details  are  referred  to  (Siriaksorn 
et  al„  2017). 
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Figure  2.  Nodal  smoothing  domain  for  the  modified  stabilized  nonconforming  nodal  integration 


4.2.  Soil  constitutive  models 

The  effective  stress  a  described  in  Section  4.1.3  is  decoupled  from  the  fluid  phase.  It  can  be 
computed  using  solid  constitutive  models  independent  of  the  fluid  phase.  Several  features  are 
necessary  to  accurately  capture  the  soil  behavior  in  the  deformation  under  penetration.  The  ability 
to  capture  large  deformation  hardening  and  softening  is  paramount.  The  ability  to  capture  the 
transition  from  dilation  to  compaction,  tensile  behavior,  rate  dependence,  and  triaxialiaty  of  the 
stress  state  can  also  be  important.  In  this  project,  two  constitutive  models  are  developed  and 
implemented  to  represent  soil  behaviors. 

4.2.1.  Non-associated  Drucker-Prager  with  damage  model 
The  Drucker-Prager  yield  function  is  described  as 

Fdp=JlT2+J3P-A  (50) 

where  J2  is  the  second  invariant  of  the  deviatoric  part  of  o ;  p  is  one-third  of  the  trace  of  a  ;  J3 

and  A  are  material  parameters;  6  is  the  effective  stress  before  degraded  by  material  damage.  The 
Drucker-Prager  material  parameters  J3  and  A  are  related  to  cohesion  c  and  friction  angle  <j>  of 
Mohr-Coulomb  by 


P  = 


2\[b  sin  (j> 


A 


3  -  sin  (f> 

2\f6c  cos  q 
3  -  sin  (j) 


(51) 

(52) 


Considering  a  non-associated  flow  rule,  in  this  project  the  plastic  potential  function  takes  the  following 
form 

G  =  y[2J\+bp-A  (53) 


where  b  can  be  related  to  the  dilatancy  angle  (//  as 
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(54) 


2>/6  sin  y/ 

b  = - 

3  -  sin  y/ 

A  damage  parameter  d  is  introduced  to  degrade  the  deviatoric  part  and  the  tensile  volumetric  part 
of  the  effective  stress  c .  It  yields  the  total  (damaged)  effective  stress  o  as 


a  =  adev  (l-d)  +  [(l-d)tr6+ +  tra  )l  (55) 

where  superscripts  dev,  +,  and  -  indicate  deviatoric,  tensile,  and  compressive  parts  of  the 
corresponding  terms,  respectively.  A  linear  evolution  of  the  damage  parameter  d  is  assumed: 


Ci(/7-c2) 

_7  ’ 


7I-C2 


(56) 


where  rj  is  the  norm  of  the  deviatoric  strain  (i.e.,  r]  =  v £dev :  £dev )  used  as  a  means  to  measure 
material  damage.  The  parameter  C2  specifies  the  initiation  point,  when  material  starts  to  degrade 
( d  =  0 ).  The  parameter  ci  specifies  the  critical  point,  when  material  is  fully  damaged  ( d  =  1 ). 


4.2.2.  Three-invariant  Soil  Constitutive  Model 

To  capture  the  transition  from  dilation  to  compaction,  tensile  behavior,  rate  dependence,  and 
triaxialiaty  of  the  stress  state,  a  three-invariant  continuum  cap  plasticity  model  is  developed  for 
large  deformation  of  soil  in  penetration  problems.  This  model  is  capable  of  simulating  multiple 
failure  mechanisms  including  loss  of  strength  under  tension  and  combined  shear  and  compaction 
yielding.  It  features  a  nonlinear  pressure-dependent  shear  yield  surface  as  well  as  a  cap  surface 
governing  inelastic  compaction  hardening  which  is  formulated  as  an  explicit  expression  of  porosity. 
It  also  accounts  for  the  differences  in  triaxial  extension  and  compression  strength.  Adapted  from 
the  previously  developed  Sandia  GeoModel  (Fossum  and  Brannon,  2004b),  the  shear  and  tension 
surfaces  are  simplified  using  a  hyperbolic  function  which  promises  a  less  computationally 
expensive  constitutive  model. 

Dynamic  problems  and  systems  exposed  to  high  strain  rates,  e.g.  soil  penetration  problems,  are 
examples  where  addressing  the  rate-dependency  factor  plays  a  key  role  in  analyzing  and  predicting 
the  system  response.  Here,  viscoplastic  regularization  by  using  an  overstress  model  of  the  Duvaut- 
Lions  type  is  adopted  in  order  to  capture  such  rate  dependency  of  the  soil  constitutive  behavior 
(Duvaut  and  Lions,  1972).  Yiscoplasticity  also  serves  to  regularize  the  constitutive  model  (Simo 
and  Hughes,  1998).  Therefore,  it  can  be  used  as  a  remedy  for  loss  of  ellipticity  issues  arising  from 
softening  or  nonassociativity. 

4.2.2. 1.  Yield  surface 

An  elliptical  tension  cap  to  better  capture  the  tensile  behavior  was  added  to  the  yield  surface  in 
(Motamedi  and  Foster,  2015).  In  this  work  modeling,  the  shear  and  tension  surfaces  are  simplified 
using  a  hyperbolic  formulation  adapted  from  (Carol  et  ah,  1997).  While  they  used  the  hyperbolic 
surface  for  yielding  along  a  discrete  surface,  the  same  concepts  can  be  adjusted  for  use  in  a 
continuum  yield  surface. 
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The  shear  yield  surface  is  given  as 


(57) 


Ff  =  yj(C-MIi)2 -(C-MZ)2 

where  /,  is  the  first  invariant  of  o ,  E  measures  the  maximum  elastic  value  of  /, ,  C  is  a  shear 

strengthen  face  parameter,  and  M  is  the  tangential  slope  of  the  shear  yield  surface  in  meridional 
stress  space  (see  Figure  3).  M  is  treated  as  constant,  but  the  other  two  parameters  may  vary. 


Figure  3.  Shear  failure  surface  Ff 


This  yield  surface  is  modified  by  an  elliptical  cap  function  Fc 


Fc{h)  =  '-H{K- 


h) 


KX-k) 


where 


(58) 


X(tc)  =  tc-RFf(tc )  (59) 

is  the  value  of  /,  where  hydrostatic  compression  occurs,  H  (.t)  is  the  Heaviside  function,  and  R 
is  a  material  parameter  that  controls  the  aspect  ratio  of  the  cap. 

The  entire  yield  surface  may  be  written  as 

f=T(p)Fl~FF,  («0) 

where  T  is  a  modifying  function  which  is  introduced  to  capture  the  difference  between  triaxial 
extension  and  compression  strength  (Fossum  and  Brannon,  2004b;  Foster  et  al.,  2005).  Thus  far, 
we  have  implemented  the  Gudehus  form  of  the  modifying  function  (Gudehus,  1973): 
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(61) 


r(A) 


2 

2 


1  1 

l  +  sin3/?H — (l-sin3/?) 

\  V  ) 

\  3V3 >3  1  f  ,  3V3/.  ^ 

v  2(72)3/2  +2(72)3/2y 


3 


y 


where  1//  is  the  constant  ratio  triaxial  extension  and  compression  strength. 


4. 2. 2. 2.  Evolution  of  plastic  variables 
(a)  Shear  and  tensile  parameters 

The  evolution  of  C  and  S  follow  a  hardening  and  softening  function.  These  strength  parameters 
may  evolve  differently  depending  on  whether  the  loading  is  compressive  or  tensile,  and  how  much 
shear  is  involved.  Hence  we  introduce  a  parameter  y  to  describe  the  evolution 


r=at(dvp}+ac(-dvp)+a 


(62) 


y  =iydt 


(63) 


where  at,  ac,and  as  are  material  constants.  The  terms  df  and  jep |  represent  the  volumetric  part 

of  the  plastic  strain  rate  and  norm  of  the  deviatoric  part  of  the  plastic  strain  rate,  respectively,  and 
(•}  are  Macaulay  brackets.  Hence,  y  is  related  to  the  plastic  multiplier  X ,  but  weighted  for  the 

type  of  deformation.  For  the  evolution  of  the  variables  C  and  S  ,  we  propose  a  hardening 
evolution  following  a  normal  function 


c(rMc,,„-c„,Kw™>:  +c„,  (64) 

s(r)=(s„ +S„,  (65) 


Here,  the  maximum  value  of  the  parameter  is  Cmax  or  Emax ,  and  residual  values  Cres  or  Eres  are 

proposed.  These  values  are  allowed  due  to  small  amounts  of  shear  and  tensile  strength  due  to 
interlocking  effects  even  in  granular  materials,  although  the  values  may  be  set  to  zero.  The  material 
constants  yc  and  yE  indicate  the  values  of  y  where  peak  strength  is  reached,  and  bc  and  bE  are 

related  to  the  initial  values  of  the  parameters  C0  and  H0  by  the  formula 


bc 

b~ 


fc  -C  ) 

max _ res 

C  -C 

V  ^0  ^ res  y 

^  max  ^ res 

V  ^0  1“l res  J 


(66) 

(67) 
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The  parameters  can  be  adjusted  so  that  this  function  predicts  hardening  and  softening,  ductile 
softening,  or  brittle  softening  depending  on  the  parameter  choices.  These  functions  are  shown  in 
Figure  4. 


C 


(a) 


(b) 


Figure  4.  Evolution  of  C  with  respect  to  y ,  (b)  Evolution  of  S  with  respect  to  y . 


(b)  Cap  parameter 

The  cap  parameter  k  is  assumed,  as  in  earlier  work,  to  be  a  function  only  of  the  volumetric  strain. 
The  evolution  of  this  variable  is  related  to  pore  collapse,  grain  crushing,  and  other 
micromechanical  phenomena.  In  order  to  couple  this  model  in  a  poromechanical  framework, 
where  permeability  changes  with  porosity,  we  need  to  track  the  changing  porosity  of  the  soil. 
Therefore,  we  reformulate  the  plastic  volumetric  deformation  as  an  explicit  expression  of  the 
porosity. 

In  small  strain,  it  can  be  shown  that  the  void  ratio  n  is  related  to  volumetric  strain  Ov  by  the 
formula 


O,  = 


n  —  n, 


o 


1  -  n 


(68) 


where  ft0  is  the  initial  void  ratio.  Similarly,  the  unloaded  porosity  (the  recovered  porosity  of  the 
material  after  elastic  unloading)  n  is  related  to  the  plastic  volumetric  strain  as 


dP  =  n  -«o 


1  —  n 


(69) 


In  large  deformation,  with  6  as  the  Eulerian  logarithmic  (Hencky)  strain,  the  formula  becomes 

(70) 


Of  =  In 


vl  ~*j 
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We  propose  a  hydrostatic  “crush  curve”  in  the  /j  -  n  plane  as 


n  =n*exp[^(X-X0)] 

(71) 

where  qx  is  a  material  constant.  Hence  as  /,  achieves  negative  large  values,  all  the  pore  space 

space  approaches  zero.  This  model  may  have  issues  with  conditioning  if  the  porosity  becomes  very 
small,  and  a  second  factor  similar  to  that  proposed  by  (Fossum  and  Brannon,  2004a)  may  improve 
the  transition  from  compaction  to  dilation. 

Recall  that 

X(*r)  :=k-RF*(k ) 

(72) 

Hence,  we  can  track  the  evolution  of  k  as 

3  f^- 

Dk  dX  dn  v  _  dl \ 

K~  dX  dn  ddvp  °v  “  dqp  dn  dX 

(73) 

dn  dX  dK 

It  can  easily  be  shown  that 

ddvp_  l-«o 
dn  (1  -n)2 

(74) 

SX  =?"»exp[9,(X  X7] 

(75) 

rsf;(k) 

dx  die 

(76) 

4.2.23.  Plastic  potential 

For  most  geomaterials,  a  plastic  potential  function  separate  from  yield  function  is  necessary  to 
prevent  the  overprediction  of  plastic  volumetric  strain.  On  shear  surface,  the  parameter  M  has  the 
greatest  influence  on  the  volumetric  dilation,  and  replacing  this  with  a  smaller  M  will  reduce  the 

dilation  on  the  shear  surface.  We  also  replace  C  with  Cg  defined  by 

C  (y)  =  (C  -C  )e~bc{r~rc)2  +C  (77) 

g  \/  )  \  max g  resg  J  resg  / 

with  the  modified  constants  Cmaxg ,  Cresg ,  C0g ,  and  yCg . 

The  plastic  potential  g  is  modified  from  the  yield  function 

g=r(p)4r2-jF*FZ  (78) 
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where 


f; = p9) 


Nonassociativity  on  the  cap  follows  (Regueiro  and  Foster,  2011),  where 

<81> 


4. 2. 2. 4.  Numerical  implementation 

Numerically,  the  RKPM  code  uses  a  finite  difference  scheme  for  time  integration,  evaluating  the 
variables  at  a  finite  number  of  discrete  points  in  time.  At  the  integration  point,  the  problem  is  strain 
driven.  Hence,  given  the  values  of  the  stress  and  internal  state  variables,  k  and  y  at  time  tn ,  along 

with  the  strain  Ql+1  at  the  current  time  ^n+1  ,  we  aim  to  find  the  updated  values  of  these  parameters 
at  time  tn+l .  In  plasticity  formulations,  this  is  typically  solved  using  a  return  mapping  algorithm. 

In  this  case,  the  equations  are  integrated  using  an  implicit  or  backward  Euler  approximation  which 
brings  first-order  accuracy  as  well  as  unconditional  stability  (at  the  integration  point  level).  First  a 
trial  stress  <r'r+1  is  computed  assuming  that  the  step  will  be  elastic 


<7 


tr 

n+ 1 


=  c„  +  Ce :  Ao 


(82) 


If  the  yield  surface  is  violated,  then  the  return  mapping  algorithm  is  invoked  to  return  the  stress  to 
the  yield  surface  and  account  for  its  evolution. 


Also,  we  employ  spectral  decomposition  of  the  stress,  as  used  in  (Borja  et  ah,  2003;  Foster  et  ah, 
2005;  Tamagnini  et  ah,  2002).  This  technique  helps  to  reduce  the  number  of  unknowns  by  working 
on  the  principal  stresses  rather  than  the  full  effective  stress  tensor.  Because  the  principal  directions, 
or  Eigenvectors,  of  the  stress  are  identical  in  the  trial  and  final  states  for  isotropic  plasticity  models, 
we  can  determine  the  directions  from  the  trial  state  and  only  solve  for  the  principal  stresses,  or 
Eigenvalues.  Adding  the  other  independent  unknowns,  the  internal  state  variables  and  the  plastic 
multiplier,  the  vector  of  unknowns  that  should  be  solved  for  at  time  tn+l  takes  the  form 


X  =  {<’"  0-“’"  Ak  Ay  AA\J 

where  cr™11  are  plastic  correctors  for  the  stress 


^corr 


=  o\ 


(83) 


(84) 
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We  solve  this  system  by  using  a  standard  Newton-Raphson  algorithm.  When  adopting  an  N-R 
algorithm,  the  iterative  calculation  results 


\rk+ 1  -yk 

A  n+l  ~  A  n+ 1  “ 


DR 

Id^A+i 


i-i 


n+1 


(85) 


where  k  + 1  denotes  the  current  iteration  number.  The  residual  vector  for  this  problem  takes  the 
form 


*(*)=■ 


AA  a. 


l  A 


AA  a 


2  A 


AA  al 


Vd(JAj 
\d(JAj 

rdg_  A 


+  cr , 


+  cri, 


+  cr 


m 


3  A 

AAhK  -Ak 

(atdvp  +ac-dvp  +asep)~  Ay 

f 


=  0 


(86) 


where  g  denotes  the  plastic  potential  function.  The  subscript  n  +  l  is  omitted  in  the  above 
equations  to  simplify  notation.  The  summation  convention  is  also  used  in  the  first  three  equations. 
The  tensor  ae  is  the  elasticity  tensor  projected  to  principal  stress  space  and  has  the  form 


a  = 


A-v'Afj, 

A 

A 


A  A 

A  +  2,/u  A 

A  A-^-A/j, 


(87) 


where  A  and  ju  are  Lame’s  first  and  second  constants,  respectively. 


The  return  mapping  algorithm  then  takes  the  form  as  outlined  in  Box  4.1. 
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Box  4.1:  Summary  of  the  return  mapping  algorithm. 


Step  1. 

Step  2. 
Step  3. 

Step  4. 
Step  5. 


Compute  <r‘r+1  =  an  +  Cc :  Ao. 

3 

Spectrally  decompose  a*+l  =  ^cr"  m{A)  . 

A=l 

Check  yielding:  is/  >  0  ? 

If  no,  set  «rn+1  =  and  exit. 

Set  X0=0  and  iterate  following  Eq.  4.31  until  the 

relative  convergence  tolerance  is  met. 

Update: 


A=1 


Kn+ 1  =  Kn  +  A** 


rn+ 1  =r„  +Ar 

and  exit. 


4. 2. 2. 5.  Parameter  fitting 

A  fitting  procedure  for  the  three-invariant  model  based  on  standard  geotechnical  tests  is  given  in 
(Fossum  and  Brannon,  2004a).  The  procedure  will  be  largely  the  same,  though  the  forms  of  the 
functions  to  be  fit  are  different.  With  regards  to  the  yield  surface,  an  additional  tensile  test  will  be 
necessary  to  accurately  fit  the  parameter  S  . 

The  hardening/softening  functions  can  still  be  fit  using,  say,  triaxial  test  data,  but  there  will  be  a 
few  more  parameters  to  tune  during  the  process.  However,  some,  such  as  yc  and  yE ,  will  have 
fairly  apparent  values,  making  them  easier  to  fit. 
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5.  Results  and  Discussion 


5.1.  Verifications 

The  three-invariant  soil  constitutive  model,  discussed  in  Section  4.2,  is  first  implemented  within 
an  FEM  code.  Using  the  three-invariant  model  and  adopting  a  specific  set  of  material  constants 
and  ignoring  the  cap,  we  are  able  to  recover  a  non-associative  Drucker-Prager  plasticity  model 
with  zero  hardening/softening.  We  used  this  characteristics  of  the  model  to  verify  the  shear  surface 
implementation,  of  the  three-invariant  model  using  a  number  of  Drucker-Prager  benchmark  tests. 

The  model  is  also  implemented  within  the  mesh-free  code.  The  two  implementations,  FEM  and 
mesh-free,  are  then  compared  by  performing  a  number  of  simulations,  including  the  uniaxial  and 
triaxial  compression,  uniaxial  tension,  and  combined  compression-shear  tests,  on  an  8-node  test 
model.  For  FEM  simulations,  a  3D  solid  cubic  finite  element  with  eight  nodes  and  eight  integration 
points  is  employed.  It  should  also  be  noted  that  the  employed  FEM  framework  uses  small  strain 
formulations  while  the  mesh-free  code  can  handle  large  deformations  and  rotations. 

The  two  frameworks,  FEM  and  mesh-free,  are  also  compared  for  the  uniaxial  and  triaxial 
compression  tests  using  a  Drucker-Prager  model  with  linear  hardening/softening.  Hence,  as  for  the 
constitutive  model,  two  cases  are  considered:  a  non-associative  Drucker-Prager  plasticity  model 
with  linear  evolution  in  the  cohesive  strength  parameter,  and  the  three-invariant  cap  plasticity 
model  regularized  by  the  Duvaut-Lions  viscoplasticity  model.  The  results  are  discussed  in  the 
following  sections. 


5.1.1.  Uniaxial  compression  test 

The  first  example,  in  FEM,  consists  of  a  one-element  model  under  uniaxial  compression,  as  shown 
in  Figure  5.  The  mesh-free  simulation  is  modeled  using  8  nodes  in  order  to  best  reproduce  the 
behavior  of  the  FEM  implementation.  In  order  to  impose  the  uniaxial  compressive  strain,  a 
displacement-controlled  approach  is  carried  out  by  applying  d  =  0.02  mm  downward  prescribed 
displacement  on  the  top  side  of  the  model,  i.e.,  the  top  four  nodes. 


1  mm 

Figure  5.  Uniaxial  compression  test.  Dimensions,  loading,  and  boundary  conditions. 
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As  for  the  constitutive  model,  a  non-associative  Drucker-Prager  plasticity  model  with  linear 
evolution  in  the  cohesive  strength  parameter  is  used  first.  The  material  properties  are  shown  in 
Table  1. 


Table  1.  Material  properties  for  the  Drucker-Prager  uniaxial  compression  test. 


Parameter 

Symbol 

Value 

Young's  modulus 

E 

72  MPa 

Poisson's  ratio 

V 

0.44 

Drucker-Prager  material  constant 

A 

0.25  MPa 

Drucker-Prager  material  constant 

P 

0.05 

Drucker-Prager  material  constant 

b 

0.025 

Hardening/softening  modulus 

Ha 

25  MPa 

For  the  Drucker-Prager  criterion  employed  here,  the  yield  function  is  expressed  as 

f  =^2T2-{K-pp)  (88) 

where  J2  is  the  second  deviatoric  stress  invariant,  a ,  cohesion-like  parameter,  (3 ,  friction-like 
parameter,  and  p  denotes  the  mean  normal  stress  defined  as 

p  =  tr(o-)/3  (89) 

with  <7  being  the  Cauchy  stress  tensor  and  tr(«)  for  the  trace  operator.  A  plastic  potential  function 
is  also  used  of  the  form 


g=^lT2-{K-bp)  (90) 

where  b  is  the  dilation  constant.  The  hardening/softening  modulus  is  used  to  demonstrate  the 
evolution  in  A  as 

A  =HaX  (91) 

hence,  forming  a  linear  hardening/softening  evolution.  The  parameter  X  represents  the  value  of 
the  plastic  strain  rate. 

The  axial  stress-strain  responses  for  the  two  frameworks  are  plotted  in  Figure  6.  A  good  agreement 
is  observed  comparing  the  two  graphs. 
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Figure  6.  Drucker-Prager  and  uniaxial  compression  test.  Compressive  axial  stress  vs.  compressive 

axial  strain  response  for  FEM  and  mesh-free. 


For  this  test,  the  three-invariant  plasticity  model  is  also  adopted  where  the  material  constants  are 
shown  in  Table  2.  The  two  frameworks  give  similar  stress-strain  responses,  as  shown  in  Figure  6. 
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Table  2.  Material  properties  for  the  three-invariant  plasticity  uniaxial  compression  test. 


Symbol 

Value 

E 

9000.0  MPa 

V 

0.15 

¥ 

0.8 

R 

1.0 

Rg 

1.0 

K0 

-50.0  MPa 

M 

0.2 

Ms 

0.2 

at 

100.0 

ac 

100.0 

as 

50.0 

^ max 

5.5  MPa 

C 

maxg 

5.5  MPa 

M 

l“l  max 

11.0  MPa 

Cres 

2.0  MPa 

c 

res  g 

2.0  MPa 

M 

l“l res 

7.0  MPa 

Yc 

1.0 

Y cg 

1.0 

Ye 

1.0 

C0 

4.5  MPa 

Cog 

4.5  MPa 

^0 

10.0  MPa 

«0 

1.0e-3 

* 

no 

1.0e-3 

E 

1.0 

T 

0.01 

26 


Figure  7.  Three-invariant  plasticity  and  uniaxial  compression  test.  Compressive  axial  stress  vs. 
compressive  axial  strain  for  FEM  and  mesh-free. 


5.1.2.  Triaxial  compression  test 

The  second  example,  in  FEM,  consists  of  a  one-element  model  under  triaxial  compression,  as 
shown  in  Figure  7.  Again,  the  mesh-free  simulation  is  modeled  using  8  nodes.  The  triaxial  test,  in 
general,  is  one  of  the  most  common  and  widely  performed  laboratory  experiments,  allowing  to 
measure  the  mechanical  properties  of  soil,  rock,  and  other  granular  materials  for  use  in  engineering 
design.  These  material  properties  include  the  angle  of  shearing  resistance,  apparent  cohesion,  and 
dilatancy  angle  among  others.  The  material  properties  for  the  Drucker-Prager  test  are  shown  in 
Table  3  and  the  parameters  have  the  same  meanings  as  in  Section  5. 1.2.1. 


1  mm 


Figure  8.  Triaxial  compression  test.  Dimensions,  loading,  and  boundary  conditions. 
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Table  3.  Material  properties  for  the  Drucker-Prager  triaxial  compression  test. 


Parameter 

Symbol 

Value 

Young's  modulus 

E 

9000.0  MPa 

Poisson's  ratio 

V 

0.15 

Drucker-Prager  material  constant 

A 

8.034  MPa 

Drucker-Prager  material  constant 

P 

0.633 

Drucker-Prager  material  constant 

b 

0.3165 

Hardening/softening  modulus 

Ha 

-100.0  MPa 

In  order  to  model  a  triaxial  compression  test,  a  given  pressure  should  be  applied  to  all  sides  of  the 
model.  A  continued  displacement  or  load  is  then  applied  in  one  direction,  maintaining  the  pressure 
on  the  other  sides.  Here,  imposing  the  triaxial  compressive  strain  in  conducted  by  using  a 
displacement-controlled  approach  where  a  prescribed  displacement  is  applied  on  the  sides  of  the 
element,  as  shown  in  Figure  7. 

The  model  is  initially  loaded  at  a  strain  rate  of  0.0004  per  second  on  all  sides  for  one  second.  This 
applies  a  confining  stress  on  the  model  similar  to  what  is  done  by  pressurizing  the  cell  fluid 
surrounding  the  specimen  in  a  laboratory  experiment.  The  strain  rate  is  then  instantaneously 
dropped  to  0  and  held  constant  for  another  second.  Then,  the  rate  is  raised  to  0.0008  per  second 
only  in  one  direction  (imposed  on  the  top  side)  for  the  rest  of  the  simulation.  Strain  rate  change 
through  the  5-second  simulation  and  the  associated  imposed  displacement,  d  ,  at  each  time  step 
are  plotted  in  Figure  8a  and  Figure  8b,  respectively.  The  axial  stress-strain  responses  are  depicted 
in  Figure  9  where  the  graphs  are  fairly  analogous. 


0.0008 


0.0004 


t  (sec) 


d  (mm) 


(a) 


(b) 


Figure  9.  Drucker-Prager  and  triaxial  compression  test,  (a)  Strain  rate  vs.  time,  (b)  Imposed 

displacement  vs.  step  number. 
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Figure  10.  Drucker-Prager  and  triaxial  compression  test.  Compressive  axial  stress  vs.  compressive 

axial  strain. 


The  three-invariant  plasticity  model  is  also  adopted  where  the  material  constants  are  shown  in 
Table  4.  For  this  model,  the  loading  conditions  are  modified  as  an  imposed  displacement  of 
d  =0.0015  mm  on  all  sides  for  the  first  second  following  a  relaxation  step  for  another  second. 
Then,  a  displacement  of  d  =  0.02  mm  is  applied  only  in  one  direction  (imposed  on  the  top  side) 
for  the  rest  of  the  simulation.  The  total  simulation  time  is  3  seconds.  The  axial  stress-strain 
responses  are  plotted  in  Figure  10.  A  similar  response  is  observed. 
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Table  4.  Material  properties  for  the  three-invariant  plasticity  triaxial  compression  test. 


Symbol 

Value 

E 

9000.0  MPa 

V 

0.15 

¥ 

0.8 

R 

2.0 

Rg 

2.0 

*0 

-50.0  MPa 

M 

0.2 

Ms 

0.2 

a, 

100.0 

ac 

100.0 

as 

50.0 

C 

max 

5.5  MPa 

cmaxg 

5.5  MPa 

M 

1“l  max 

11.0  MPa 

Cres 

2.0  MPa 

c 

res  g 

2.0  MPa 

1—1 res 

7.0  MPa 

Yc 

1.0 

Y cg 

1.0 

Ye 

1.0 

Co 

4.5  MPa 

Co, 

4.5  MPa 

M 

^0 

10.0  MPa 

n0 

1.0e-l 

* 

n0 

1.0e-l 

<h 

10.0 

T 

0.01 
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Figure  11.  Three-invariant  plasticity  and  triaxial  compression  test.  Compressive  axial  stress  vs. 

compressive  axial  strain. 


5.1.3.  Uniaxial  extension  test 

The  next  two  simulations  are  only  performed  using  the  three-invariant  model.  This  example,  in 
FEM,  consists  of  a  one-element  model  under  uniaxial  extension,  as  shown  in  Figure  1 1 .  The  mesh- 
free  simulation  is  modeled  using  8  nodes  as  before.  By  performing  this  test,  we  aim  to  check  the 
tension  cap  of  the  model.  In  order  to  impose  the  uniaxial  tensile  strain,  a  displacement-controlled 
approach  is  carried  out  by  applying  d  =  0.01  mm  upward  prescribed  displacement  on  the  top  side 
of  the  model,  i.e.,  the  top  four  nodes. 


1  mm 

Figure  12.  Uniaxial  extension  test.  Dimensions,  loading,  and  boundary  conditions. 
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The  material  constants  for  the  three-invariant  model  are  identical  to  those  in  Table  2,  for  uniaxial 
compression  test,  except  for  M  =  0.05  which  is  chosen  in  order  to  test  nonassociativity.  Figure 

12  compares  the  axial  stress-strain  responses  for  the  two  frameworks.  The  graphs  are  comparable. 


Figure  13.  Three-invariant  plasticity  and  uniaxial  extension  test.  Axial  stress  vs.  axial  strain  for 

FEM  and  mesh-free. 


5.1.4.  Combined  compression-shear  test 

The  last  simulation  includes  a  compression  step  followed  by  a  shear  loading  mainly  to  verify  the 
model  under  a  rotation  of  the  principal  stress  axes.  Figures  13a  and  Figure  13b  provide  the  loading 
and  boundary  conditions  of  the  two  steps.  During  compression,  a  displacement-controlled 
approach  is  carried  out  by  applying  dc  =  0.01  mm  downward  prescribed  displacement  on  the  top 
side  of  the  model.  Next,  the  top  face  is  constrained  to  move  in  vertical  direction,  and  the  top  four 
nodes  are  moved  horizontally  by  applying  ds  -  0.02  mm. 
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(a)  (b) 

Figure  14.  Combined  compression-shear  test.  Dimensions,  loading,  and  boundary  conditions  for  (a) 

compression  step,  (b)  shear  step. 


The  material  constants  for  the  three-invariant  model  are  identical  to  those  in  Table  4,  for  triaxial 
compression  test,  except  for  qx  =  0.001 .  The  shear  stress-strain  responses  for  the  two  frameworks 
are  depicted  in  Figure  14.  The  graphs  are  comparable. 


Figure  15.  Three-invariant  plasticity  and  combined  compression-shear  test.  Shear  stress  vs.  shear 

strain  for  FEM  and  mesh-free. 


5.1.5.  Discussion 

Throughout  Section  5.1.2,  given  the  formulations  developed  in  Section  4.2,  most  of  the  salient 
characteristics  of  the  proposed  smooth  three-invariant  cap  model  have  been  tested  and  verified  by 
a  number  of  numerical  examples. 

The  above  examples  demonstrate  some  of  the  important  features  of  the  new  model  in  numerical 
frameworks  as  well  as  verify  the  two  implementations.  With  regards  to  the  latter,  the  model  has 


33 


been  first  verified  against  the  well-known  Drucker-Prager  model  to  the  extent  that  the  models  are 
comparable.  The  two  models  matched  up  to  machine  round  off  error.  Comparing  the  finite  element 
and  meshfree  implementations  of  the  new  model,  one  can  see  the  two  implementations  are  nearly 
identical.  Small  differences  are  due  to  the  dynamic  implementation,  albeit  with  small  mass,  and 
explicit  global  time  stepping  of  the  meshfree  code,  as  well  as  a  different  number  of  time  steps.  In 
addition,  the  FEM  formulation  uses  a  small  strain  assumption  while  the  meshfree  accounts  for 
large  strains. 

The  examples  also  verify  that  the  intended  features  of  the  models  work  under  uniaxial  and 
multiaxial  loading  conditions.  The  hardening  and  softening  behavior  of  the  model  is  evident  in 
the  examples,  as  is  the  rate  dependence.  The  new  tension  cap  is  also  performing  as  expected. 
Though  not  discussed  here,  the  compression  cap  and  the  triaxiality  functions  have  also  been 
tested  in  the  finite  element  framework. 


5.2.  Ad  hoc  projectile  penetration  in  soil  simulations 

To  demonstrate  the  capability  of  the  developed  two-field  meshfree  code,  ad  hoc  projectile- soil 
penetration  tests  with  different  impact  angles  are  set  up  as  described  in  Table  5.  The  projectile’s 
geometry  and  soil’s  dimensions  are  set  up  to  be  the  same  and  an  initial  velocity  of  150  m/s  for 
the  projectile  is  considered  for  all  cases.  The  dimensions  of  the  soil  domain  are  16mxl2m  (the 
impact  surface)  with  depth  of  8m.  They  are  chosen  large  enough  to  ensure  that  reflected  waves 
from  the  boundaries  have  minimal  influences  on  the  penetration  area.  To  mimic  a  semi-infinite 
soil  domain,  fixed  boundary  conditions  are  employed  on  all  surfaces  of  the  soil  domain,  except 
the  traction-free  surface  (top  impact  surface). 
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Table  5.  Problem  setup  of  projectile  penetration  in  soil 


Impact  Angle 

(Degrees) 


90° 


80° 


70° 


60° 


45° 


30° 


A  cylindrical  projectile  with  a  radius  of  0.2m  and  a  height  of  2.5m  is  selected  to  mimic  the 
training  munitions  dropped  on  the  test  site.  The  detailed  projectile’s  dimensions  are  shown  in 
Figure  15. 
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Figure  16.  Dimension  of  projectile 

Although  the  mesh  is  not  required  for  the  meshfree  simulation,  a  standard  FEM  mesh  software  is 
used  to  create  the  point  discretization  for  the  projectile  and  the  soil  domain.  The  mesh 
configurations  for  the  projectile  and  the  soil  are  shown  in  is  selected  as  shown  in  Figure  16.  The 
nodal  distance  is  fairly  uniform  for  the  whole  model  (projectile  and  soil),  0.151m  on  average. 
This  spacing  keeps  a  good  balance  between  accuracy  and  computational  efficiency. 


Figure  17.  Discretization  for  the  projectile  and  the  soil 


In  this  simulation,  the  projectile  is  modeled  as  linear  elastic  with  properties  given  in  Table  6.  The 
soil  is  modeled  with  the  non-associative  Drucker-Prager  model  with  damage  (Section  4.3)  under 
the  two-field  semi-Lagrangian  RK  formulation  (Section  4.2).  The  material  properties  of  the  soil 
model  are  given  in  Table  7. 
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Table  6.  Material  properties  of  projectile 


Linear  Elastic  Material 

Young’s  modulus  (E) 

2E+H(£) 

Poisson’s  ratio  (v) 

0.25 

Density  (p) 

8050  (4) 

Mass  proportional  damping 

0.05 

Table  7.  Material  properties  of  soil 


Drucker-Prager  geomaterial  properties 

Young’s  modulus  (E) 

2E  +  8  (4) 

Poisson’s  ratio  (v) 

0.2 

Friction  (/?) 

0.16 

Hardening 

IE +  8  (4) 

Cohesion  strength  (a) 

2E  +  4  (4) 

First  parameter  for  damage  accumulation  function 

0.05 

Second  parameter  for  damage  accumulation  function 

1.0 

Density  (p) 

2000  (4) 

m6 

Mass  proportional  damping 

0.05 

The  deformation  snapshots  of  the  penetration  for  each  case  are  given  in  Figures  18,  20,  22,  24, 
26,  28.  The  color  in  the  plots  indicates  the  level  of  damage.  The  time  history  data  for  the 
impact/penetration  process  in  the  first  0.2  second  are  reported  in  Figures  19,  21,  23,  25,  27,  29. 

In  each  case,  the  results  are  presented  using  five  different  plots: 

(a)  Trajectory: 

This  plot  shows  the  projectile-tip’s  node  path  during  impact  and  penetration.  The 
maximum  depth  of  penetration  can  be  observed. 

(b)  Horizontal  displacement  time  history: 

This  plots  shows  the  horizontal  position  of  projectile’s  tip  with  respect  to  time. 

(c)  Vertical  displacement  time  history: 

This  plots  shows  the  vertical  position  of  projectile’s  tip  at  each  time  step.  On  this  plot,  the 
maximum  penetration  depth  and  the  time  it  occurs  are  shown. 

(d)  Horizontal  velocity  time  history: 

This  plot  shows  the  horizontal  nodal  velocity  at  different  times. 

(e)  Vertical  velocity  time  history: 

This  plot  shows  the  vertical  nodal  velocity  at  different  times. 
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Table  8  shows  the  maximum  penetration  depth  for  all  cases: 


Table  8.  Projectile’s  maximum  penetration  depth  for  different  impact  angles 


Case  No.  (Impact 
Angle) 

Maximum  Penetration  Depth 

(m) 

#1  (90°) 

3.00 

#2  (80°) 

2.96 

#3  (70°) 

2.40 

#4  (60°) 

2.2 

#5  (45°) 

1.08 

#6  (30°) 

0.08 

Figure  18.  Deformation  of  penetration  in  soil  (impact  angle  90°) 
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0.5 


P  rojecti  le  Ti  p  -N  od  e  T  raj  ecto  ry 


Tip  Node  X  position  vs.  Time  Tip  Node  Y  position  vs.  Time 


Figure  19.  Trajectory  and  time  history  of  penetration  (impact  angle  90°).  (a)  tip  trajectory,  (b) 
horizontal  displacement  history,  (c)  vertical  displacement  history,  (d)  horizontal  velocity  history,  (c) 

vertical  velocity  history 
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Figure  20.  Deformation  of  penetration  in  soil  (impact  angle  80°) 
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P  rojecti  le  Ti  p  -N  od  e  T  raj  ecto  ry 


Tip  Node  X  position  vs.  Time  Tip  Node  Y  position  vs.  Time 


Figure  21.  Trajectory  and  time  history  of  penetration  (impact  angle  80°).  (a)  tip  trajectory,  (b) 
horizontal  displacement  history,  (c)  vertical  displacement  history,  (d)  horizontal  velocity  history,  (c) 

vertical  velocity  history 
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Figure  22.  Deformation  of  penetration  in  soil  (impact  angle  70°) 
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P  rojecti  le  Ti  p  -N  od  e  T  raj  ecto  ry 


Tip  Node  X  position  vs.  Time  Tip  Node  Y  position  vs.  Time 


Figure  23.  Trajectory  and  time  history  of  penetration  (impact  angle  70°).  (a)  tip  trajectory,  (b) 
horizontal  displacement  history,  (c)  vertical  displacement  history,  (d)  horizontal  velocity  history,  (c) 

vertical  velocity  history 
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Figure  24.  Deformation  of  penetration  in  soil  (impact  angle  60°) 


45 


Projectile  Tip-Node  Trajectory 


Tip  Node  X  position  vs.  Time  Tip  Node  Y  position  vs.  Time 


Figure  25.  Trajectory  and  time  history  of  penetration  (impact  angle  60°).  (a)  tip  trajectory,  (b) 
horizontal  displacement  history,  (c)  vertical  displacement  history,  (d)  horizontal  velocity  history,  (c) 

vertical  velocity  history 
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Figure  26.  Deformation  of  penetration  in  soil  (impact  angle  45°) 
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P  rojecti  le  Ti  p  -N  od  e  T  raj  ecto  ry 


Tip  Node  X  position  vs.  Time  Tip  Node  Y  position  vs.  Time 


Figure  27.  Trajectory  and  time  history  of  penetration  (impact  angle  45°).  (a)  tip  trajectory,  (b) 
horizontal  displacement  history,  (c)  vertical  displacement  history,  (d)  horizontal  velocity  history,  (c) 

vertical  velocity  history 
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Figure  28.  Deformation  of  penetration  in  soil  (impact  angle  30°) 


49 


Tip  Node  Y  position  vs.  Time 


Tip  Node  X  position  vs.  Time  Projectile  Tip-Node  T rajectory 


Figure  29.  Trajectory  and  time  history  of  penetration  (impact  angle  30°).  (a)  tip  trajectory,  (b) 
horizontal  displacement  history,  (c)  vertical  displacement  history,  (d)  horizontal  velocity  history,  (c) 

vertical  velocity  history 
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5.2.1.  Conclusions 


Since  experimental  data  for  validations  of  penetration  into  soil  is  limited,  a  complete  set  of  data 
with  different  parameters  have  not  found  in  the  public  literature.  The  ad  hoc  simulations  in  this 
section  is  aimed  to  test  the  performance  of  the  newly  developed  and  implemented  two-field 
semi-Lagrangian  RK  code. 

Despite  no  corresponding  experimental  data  available  for  validation,  the  simulations 
demonstrate  firstly  that  stable  numerical  solutions  can  be  obtained.  The  stability  is  a  particular 
concern  in  this  type  of  simulations.  Instability  may  be  caused  by  nodal  domain  integration,  the 
approximations  for  the  displacement  and  pressure  fields,  explicit  temporal  integration,  and/or 
combination  of  abovementioned.  The  current  implementation,  with  the  MSNNI,  linear  RK  for 
both  displacement  and  pressure  field,  and  explicit  temporal  integration,  provides  an  effective 
means  to  simulate  the  penetration  problems. 

Secondly,  the  basic  inelastic  impact  mechanisms  are  captured  in  the  simulations:  the  maximum 
penetration  depth  reduces  as  the  impact  angle  reduces,  and  reflection  angle  is  smaller  than  the 
impact  angle.  The  soil  deformation  during  the  penetration  also  shows  material  splashing  on  the 
free  surface,  which  is  commonly  seen  in  similar  experiments. 

One  interesting  observation  in  the  numerical  test  is  that  the  projectile  is  bouncing  back  even 
if  the  impact  angle  is  close  to  90°.  This  may  be  attributed  to  finite  boundary  conditions,  over 
stiffness,  low  cohesion,  and  damping  in  soil  properties.  On  the  other  hand,  the  observation  may 
indeed  reflect  the  physical  penetration  process  happened  in  the  experiments.  In  such  case,  the 
final  penetration  depth,  after  splashing  material  and  projectile  settle  down,  is  anticipated  to  be  less 
than  the  maximum  penetration  depth.  More  studies  and  designed  experiments  are  needed  to 
validate  the  numerical  simulations. 


5.3.  Simulations  of  spherical  ball  drop  test 

The  motivation  behind  this  simulation  is  constructing  a  numerical  model,  by  using  the  calibrated 
material  properties  and  the  Drucker-Prager  constitutive  model  used  in  the  previous  section,  to 
reproduce  and  validate  the  experimental  results  reported  by  (Seguin  et  al.,  2008).  For  this 
purpose,  a  steel-made  sphere  ball  falling  under  its  own  weight  on  an  unbounded  grain  medium  is 
modeled  and  its  maximum  penetration  has  been  evaluated.  The  term  unbounded  states  the  fact 
that  the  grain  container  dimensions  are  chosen  to  be  large  enough  to  decrease  the  boundary 
effects  on  the  results  to  its  minimum. 

The  experimental  setup  for  this  test  is  shown  in  Figure  30: 
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Figure  30.  Experimental  setup  (Seguin,  et  al  2008) 


Variables  shown  in  the  Figure  30  descried  as  below: 

S\  Maximum  penetration  depth. 

h:  Initial  distance  between  the  sphere  ball’s  lowest  point  and  grain  container’s  surface. 
H:  Total  drop  height  after  penetration. 
d :  Sphere  ball  diameter  (19  (mm)  for  this  test) 

D:  cylindrical  container’s  diameter  (190  (mm)  for  this  test) 
b:  cylindrical  container’s  height  (150  (mm)  for  this  test) 


In  this  simulation,  three  different  initial  heights  (h)  of  100, 300,  and  500  (mm)  are  tried.  For 
computational  efficiency,  for  all  three  cases,  the  ball  is  located  at  the  distance  of  50  (mm)  and 
corresponding  initial  velocities  are  assigned  to  the  ball  nodes  for  each  case.  For  this  purpose,  the 
difference  of  actual  total  height  (h)  and  the  assumed  initial  numerical  height  are  used  by  the 
formula  below, 


v  =  y[2g(h  —  0.05) 


where  g  is  the  ground  acceleration  of  9.806 


jj)  and  v  is  the  ball’s  initial  vertical  speed. 
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Table  9  shows  the  initial  velocity  for  three  case.  The  granular  medium  consists  of  glass  beads 
with  properties  presented  in  Table  10  and  the  steel-made  ball  linear  elastic  properties  are  as 
given  in  Table  11. 


Table  9.  Initial  position  &  Velocity  for  spherical  ball 


Drop  Height 
(mm) 

Height  Difference 
(mm) 

Initial  Velocity 

(mm/s) 

500 

450 

2970.85 

300 

250 

2214.34 

100 

50 

990.28 

Table  10.  Drucker-Prager  parameters  &  properties 


Drucker-Prager  granular  glass  beads  properties 

Young’s  modulus  (£) 

2.0486 E  +  8  (-^) 

Poisson’s  ratio  (v) 

0.3227 

Friction  (/?) 

0.1003 

Hardening 

0.0 

Cohesion  strength  (a) 

0.0 

First  parameter  for  damage  accumulation  function 

0.05 

Second  parameter  for  damage  accumulation  function 

1.0 

Density  (p) 

2500  (%) 

Mass  proportional  damping 

0.05 

Table  11.  Spherical  ball  mechanical  properties 

Linear  Elastic  Material  (steel  ball) 

Young’s  modulus  (E) 

2E+11(£) 

Poisson’s  ratio  (v) 

0.25 

Density  (p) 

7800  (^r) 

Mass  proportional  damping 

0.05 
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An  average  mesh  size  of  3.5  (mm)  of  hexahedral  element  is  used  for  both  ball  and  grain  parts. 
Here  the  mesh  configurations  are  shown  for  both  parts.  Figure  31  shows  the  mesh  configuration 
for  this  model. 


Figure  31.  (a)  Spherical  ball  mesh  configuration,  (b)  glass  beads  medium  mesh  configuration, 
(c)  Initial  ball’s  position  in  the  numerical  model 
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Figures  32  to  34  show  different  fames  of  the  test.  Figures  on  the  right-hand  side  show 
approximately  an  equal  time  step  of  the  simulation.  As  it  is  clear,  the  higher  the  drop  height 
becomes,  the  more  severe  the  impact  is  observed. 


Figure  32.  Casel;  500(mm)  drop  height 


Figure  33.  Case2;  300(mm)  drop  height 
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Final  penetration  depth  for  all  three  cases  are  shown  in  the  Table  12: 


Table  12.  .normalized  total  drop  height  and  penetration  depth 


S  (mm) 

h  (mm) 

H  (mm) 

H/d 

S/d 

S/d 

(Experiment) 

1.1207 

100 

115.8 

6.1 

0.83 

1.4  ±0.2 

5.6147 

300 

320.4 

16.86 

1.07 

2.4  ±  0.2 

7.6883 

500 

526.3 

27.7 

1.38 

2.85  ±  0.2 

As  shown  in  the  table  above,  increasing  the  ball’s  release  height  would  result  in  a  higher 
penetration  depth.  However,  getting  better  agreement  with  the  experimental  data,  requires  a 
better  means  or  more  experimental  testing  to  calibrate  material  constants.  Figure  35  shows 
penetration  depth  along  time  for  the  three  simulated  cases. 


Penetration  Depth  vs.  Time 


Figure  35.  Penetration  history  for  drop  test 


5.3.1.  Discussion 

The  description  of  the  material  parameters  of  the  glass  beads  in  the  Seguin,  2008  paper  is  not 
complete.  We  have  used  estimates  of  the  elastic  parameters  and  friction  angle  from  the  literature. 
Considerable  variation  in  the  properties  may  exist,  especially  in  the  initial  porosity,  which 
greatly  affects  fiction  angle.  The  measured  elastic  modulus  also  varies  significantly  for  glass 
beads,  and  the  fitting  procedure  for  the  bulk  elastic  modulus  from  this  is  approximated.  Finally, 
the  Drucker-Prager  model  is  approximate  and  does  not  capture  perfectly  all  of  the  behaviors  of 
the  beads  during  impact. 

As  we  improve  the  validation,  we  will  use  a  few  simulations  to  better  fit  the  bulk  elastic  and 
plastic  parameters.  Those  parameters,  in  turn,  will  be  used  perform  other  simulations  not  used  in 
the  fitting,  to  validate  that  the  model  can  capture  physical  phenomena  with  accurate  input. 
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6.  Conclusions  and  Implications  for  Future  Research/Implementation 

6.1.  Conclusions: 

The  two-field  (displacement-pressure)  formulation  based  on  Biot  theory  has  been  developed  and 
implemented  under  the  semi-Lagrangian  RK  framework,  where  displacement  and  pressure  field  are 
independently  approximated  by  the  semi-Lagrangian  RK  shape  functions.  Numerical  schemes 
originally  designed  for  the  single-field  formulation  have  been  modified  and  implemented  for  the  two- 
field  formulation,  including  the  modified  stabilized  non-conforming  nodal  integration  for  the  domain 
integration,  stress  update,  and  kernel  contact  algorithms.  The  central  difference  and  forward  Euler 
temporal  integration  schemes  have  been  applied  to  the  displacement  and  pressure  fields,  respectively, 
in  the  two-field  formulation,  leading  to  an  explicit  time  marching  scheme. 

To  present  the  soil  behavior,  two  constitutive  models  are  updated  and  implemented  in  the  semi- 
Lagrangian  framework:  the  Drucker-Prager  plasticity  with  damage  model  and  three  invariant 
viscoplasticity  model.  In  the  former  model,  a  single  damage  parameter  is  introduced  to  degrade  the 
deviatoric  and  tensile  parts  of  the  effective  stress,  providing  a  simple  means  to  represent  material 
damage  and  softening. 

A  three-invariant  viscoplasticity  model  has  been  developed  to  effectively  integrate  tensile,  shear, 
and  compressive  behavior.  The  evolution  of  volumetric  plastic  strain  has  been  explicitly 
connected  to  the  void  ratio,  allowing  the  model  to  be  integrated  with  a  poromechanical 
framework.  Novel  hardening/softening  laws  to  have  been  added  to  characterize  strengthening 
and  weakening  in  different  loading  regimes,  and  regularized  the  softening  using  viscoplasticity. 
The  model  also  accounts  for  rate  effects,  differences  in  strength  in  triaxial  extension  and 
compression,  compression  hardening,  and  other  effects.  An  efficient  implementation  using  the 
spectral  decomposition  has  been  employed  to  reduce  the  cost  of  the  complicated  model. 

The  model  has  been  verified  against  a  Drucker-Prager  model,  and  the  meshfree  and  finite 
element  versions  of  the  model  have  been  verified  against  each  other  to  ensure  proper 
implementation.  The  behaviors  of  the  model  have  been  demonstrated  in  a  numerical  framework 
using  reasonably  simple  example  problems. 

The  developed  two-field  meshfree  code  has  been  employed  to  simulate  penetration  process  into  soil 
and  predict  the  final  penetration  depth  under  different  penetration  angles.  In  the  ad  hoc  simulations, 
the  numerical  results  show  that  the  maximum  penetration  depth  varies  from  3m  to  lm  with 
penetration  angles  ranging  from  90°  to  45°.  The  deformation  of  soil,  e.g.  soil  splashing  on  the  free 
surface,  reflects  the  experiment  observations  after  impact.  The  developed  meshfree  code  has  been 
used  to  simulate  a  spherical  ball  impact  a  granular  medium  (Seguin,  2008).  The  numerical  results, 
at  the  current  stage,  does  not  agree  with  the  experimental  data  due  to  lack  of  data  to  calibrate 
material  constant  in  the  simulation  and  the  reasons  discussed  in  the  Section  5.3.1.  Nevertheless, 
the  numerical  results  show  that  the  penetration  depth  increases  with  the  increase  of  drop  height 
and  that  the  penetration  depth  increase  from  the  drop  height  of  100mm  to  300mm  is  higher  than 
that  from  the  drop  height  of  300mm  to  500mm,  which  agree  with  experimental  observations. 

Based  on  the  numerical  studies  and  observations  in  the  preliminary  simulations,  suggestions  to 
improve  the  performance  of  the  meshfree  code  and  soil  constitutive  models  and  the  calibration  of 
material  properties  are  given  in  the  following  section. 
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6.2.  Future  Research/Implementation: 

6.2.1.  Infinite  boundary  in  meshfree  framework 

The  physical  domain  of  soil  is  best  described  as  semi-infinite.  The  computational  domain, 
however,  is  limited  to  a  finite  domain  for  efficiency.  In  the  current  study,  the  fixed  boundary 
condition  was  used  for  the  surfaces  that  are  connected  with  surrounding  soils.  The  fixed 
boundary  condition,  when  hit  by  elastic  waves,  generates  unphysical  reflection  waves  that 
disturb  the  desired  solutions  near  the  impact  areas.  To  avoid  the  reflected  wave  due  to 
boundaries,  the  computation  domain  for  soil  in  the  project  was  constructed  large  enough  to  avoid 
the  reflected  waves  interacting  with  penetration  processes.  However,  when  the  total  simulation 
time  increases,  the  reflected  waves  eventually  interfere  penetration  processes. 

A  more  efficient  way  to  deal  with  the  issues  is  implementation  of  the  “infinite  element” 
(Zienkiewics  et  al.  1983)  in  the  meshfree  formulation.  Through  the  explicit  or  implicit 
enrichment  (Chen  et  al.  2017),  the  decaying  function  (Zienkiewics  et  al.,  2013)  due  to  infinite 
domain  can  be  embedded  in  the  RK  approximation  for  handling  the  semi-infinite  characteristics 
of  the  problem.  The  infinite  element  can  be  developed  for  both  the  single-field  or  two-field 
formulations. 

6.2.2.  Gradient  enhanced  stability  and  quasi  linear  formulation 

The  nodal  integration,  one  integration  point  per  node,  is  computationally  efficient  in  the  Galerkin 
meshfree  formulation.  This  integration  scheme,  however,  has  known  accuracy  and  stability  issues. 
Although  the  Modified  Stablized  Nonconforming  Nodal  Integration  has  been  implemented  in  the 
current  meshfree  code,  numerical  instability  may  occur  under  certain  conditions.  The  numerical 
stability  issues  may  be  attributed  to  many  reasons,  such  as  the  under  integration  in  the  Galerkin 
formulation,  losing  neighbor  nodes  needed  to  form  the  linear  basis  functions,  Eulerian  kernels,  and  u- 
p  formulation.  To  address  the  issue  associated  with  under  integration,  the  advanced  stability  (Hillman 
and  Chen,  2016)  or  (Wu  et  al.,  2016)  can  be  implemented  in  the  meshfree  code  and  be  extended  to 
two-field  formulation.  The  quasi-linear  reproducing  kernel  (Yreux  and  Chen,  2017)  can  be 
implemented  in  the  two-field  formulation  to  remedy  the  issues  of  losing  neighbor  points,  and 
ultimately  to  increase  the  stability  of  the  meshfree  framework. 


6.2.3.  Regularization  for  damaged/softened  material  behaviors 

Strain  softening  or  damage  in  the  constitutive  model  causes  strain  localization  and  the  numerical 
solutions  exhibit  pathological  discretization  sensitivity.  The  numerical  solution  does  not 
converge  with  the  refinement  of  discretization  due  to  loss  of  ellipticity  in  the  boundary  value 
problem.  Many  regularization  techniques  have  been  proposed  to  remedy  the  strain  localization 
and  discretization  sensitivity  issues  (de  Borst,  2004;  Besson,  2011).  One  of  most  effective 
methods  is  to  introduce  the  strain  gradient  in  the  variational  equation.  However,  the  strain 
gradient  requires  the  second  derivatives  of  approximation,  making  computation  very  costly. 
Instead  of  the  direct  derivatives,  the  implicit  derivatives  can  be  efficiently  obtained  under  the 
graduate  RK  approximation.  Therefore,  the  implicit  strain  gradient  is  suggested  to  be  included  in 
the  current  two-field/multi-field  meshfree  code. 
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6.2.4.  Constitutive  Model  Improvements 

While  the  constitutive  model  implemented  reproduces  most  relevant  behaviors  of  the  soil  well, 
further  developments  can  improve  both  the  fidelity  of  the  model  and  performance  of  the 
numerical  algorithm.  The  following  are  the  key  improvements  of  the  model  proposed  to 
accurately  and  efficiently  reproduce  soil  behavior. 

Currently,  the  frictional  strength  of  the  material  is  considered  constant.  In  reality,  as  materials  are 
compacted,  interlocking  effects  and  frictional  strength  increase.  We  will  modify  the  parameter  M 
to  be  function  of  the  porosity,  increasing  as  the  material  is  compacted  and  decreasing  as  it 
dilates. 

Plastic  compaction  can  also  begin  in  the  shear  strengthening  regime.  In  order  to  capture  this 
behavior,  we  will  disconnect  the  plastic  potential  variable  k?  from  the  variable  k.  This  update 
will  allow  the  tracking  of  the  volumetric  plastic  strain  separately  from  the  yield  surface.  A 
similar  evolution  law  for  both  k  and  k?  will  still  be  applicable,  and  hence  this  modification  can 
be  made  with  rather  a  rather  minor  update  to  the  model. 

Finally,  we  will  incorporate  nonlinear  elasticity  into  the  constitutive  model.  Especially  clay  soils 
exhibit  a  nonlinear  elastic  response.  While  the  effects  of  nonlinear  elasticity  can  be 
approximately  captured  by  modifying  the  plasticity  model,  in  complicated  loading  scenarios  that 
involved  unloading  it  is  better  incorporate  nonlinear  elasticity.  Currently,  we  are  implementing  a 
geometrically  nonlinear  model  that  is  materially  linear  between  the  Kirchhoff  stress  and  Hencky 
strain,  following  (de  Souza  Neto,  et  al.,  2008).  Several  nonlinear  elastic  models  have  been 
proposed  for  soils.  Borja  (2013)  reviews  several  of  those  models,  and  we  will  follow  those 
models,  modifying  as  necessary. 

The  modified  constitutive  model  will  be  integrated  into  a  partially  saturated  soil  framework. 
Currently,  we  have  implemented  a  fully  saturated  soil  model  in  conjunction  with  the  plasticity 
model.  The  extension  the  partially  saturated  case  is  easier  in  some  respects,  in  that  volumetric 
locking  issues  are  not  likely  to  be  a  problem.  There  are  a  number  of  approaches  to  modeling 
partially  saturated  media,  see  discussions  in  Coussy,  et  al.,  2004  and  Gens,  et  al.  2006  for  details. 
A  solid  plasticity  model  can  be  incorporated  into  a  partially  saturated  framework  following  Borja 
and  White  (2010)  using  a  generalized  effective  stress.  This  will  be  the  initial  approach  for  the 
new  implementation. 

6.2.5.  Numerical  Implementation  Improvements  for  soil  models 

Furthermore,  we  will  investigate  methods  to  improve  the  efficiency  and  robustness  of  the 
numerical  algorithm.  Because  the  return  mapping  perpendicular  to  the  hydrostatic  axis  is  nearly 
the  same  in  the  trial  and  final  states,  it  may  be  possible  to  use  the  trial  value  of  the  stress  gradient 
of  this  part  of  the  plastic  potential,  leading  to  a  scalar  equation  to  be  solved  rather  than  a  matrix 
equation.  Furthermore,  because  this  part  of  the  return  is  bounded  radially,  it  may  not  affect 
stability.  We  will  formally  investigate  the  accuracy  and  stability  implications  of  this  approach,  as 
well  as  the  gains  in  efficiency. 

Another  alternative  is  to  implement  a  semi-implicit  (e.g.  Belytschko  et  al.,  2000)  or  explicit 
implementation,  since  the  global  time  stepping  scheme  is  already  explicit.  This  will  not  affect  the 
order  of  accuracy  but  may  somewhat  reduce  the  time  step  needed  for  stability. 
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6.2.6.  Multiscale  material  modeling  for  calibration  of  material  properties 

One  of  the  main  causes  of  discrepancy  between  the  numerical  results  and  experimental  data  is  the 
material  constant  calibration.  While  some  well-controlled  experiments  report  data  for  penetration 
process,  the  necessary  tests  for  calibration  of  material  constant  for  soil  constitutive  models  may  not 
be  provided  or  conducted.  Moreover,  in  the  meso-  or  micro-  scale,  the  continuum  theory  may  not  be 
applicable  for  granular  materials.  A  multiscale  framework  is  therefore  recommended.  In  the  meso-  or 
micro-scale,  the  discrete  element  method  (DEM)  can  be  applied  to  obtain  the  mechanical  behavior  of 
the  granular  materials  within  a  representative  volume  element  (RVE).  Then  the  stress-strain 
relationship  of  the  granular  material  occupying  the  RVE  can  be  homogenized  and  passed  to  the 
continuum  models  in  the  macroscale  (Andrade  an  Tu,  2009;  Andrade,  et  al.,  2011;  Ren  et  al.,  201 1). 


6.2.7.  Parametric  studies  and  suggestion  of  penetration  equations 

Once  the  accuracy  and  robustness  of  the  meshfree  is  verified  and  validated,  parametric  studies 
can  be  conducted  to  understand  and/or  predict  how  different  penetration  parameters  affect  the 
final  penetration  depth  for  various  soil  types.  Ultimately,  penetration  equations  can  be  suggested 
based  on  simulations. 
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